Meta-analysis of the Interoceptive Accuracy Scale (IAS)

Study 1

Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Data Preparation

  • Murphy (2020)
    • https://osf.io/3m5nh
Code
df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))
  • Gaggero (2021)
    • https://osf.io/5x9sg
Code
load("../data/Gaggero2021/data.rda")
df2 <- data |> 
  mutate(Gender = as.character(gender)) |> 
  mutate(across(starts_with("IAS "), as.numeric)) |>
  rename(
    Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
    Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
    Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
    Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
    Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)
  • Campos (2022) - Study 1
    • https://osf.io/j6ef3
Code
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )
  • Todd (2022)
    • https://osf.io/ms354/
Code
df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
  rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )
  • Arslanova (2022)
    • https://osf.io/mp3cy/
    • note: Final.xlsx includes information about which participants passed the attention checks
Code
df5_attention <- readxl::read_excel("C:/Users/asf25/Documents/studies/InteroceptionIAS/data/Arslanova2022/Participants2.xlsx", 
    sheet = "FINAL") |>
  filter(ActorHR == 1) |>
  select(Subj.ID, Gender, Age) |>
  mutate(Gender = as.numeric(Gender), Age = as.numeric(Age)) |>
  filter(!is.na(Age)) |> #error on the documentation of this participant - attention check failed but not noted
  mutate(Gender = case_when(
    Gender== 0 ~ "Male",
    Gender== 1 ~ "Female" # based on paper reporting 65 women
  ))
New names:
• `Why No?` -> `Why No?...4`
• `Why No?` -> `Why No?...6`
• `` -> `...12`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`
Code
df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx")  |>
  filter(str_detect(Question.Key, "IAC_")) |> 
  filter(str_detect(Question.Key, "-quantised")) |> 
  pivot_wider(names_from = Question.Key, values_from = Response) |>
  rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |> 
  rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |> 
  rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
         Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast, 
         Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar, 
         Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |> 
  rename(Subj.ID = "Participant.Public.ID") |>
  select(-starts_with("Participant")) 

df5 <- df5 |>
  filter(Subj.ID %in% df5_attention$Subj.ID) |>
  select(-Subj.ID) |>
  mutate(across(everything(), as.numeric))
  • Brand (2022)
    • https://osf.io/xwz6g/
Code
df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |> 
  select(-SERIAL) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  • Brand (2023)
    • https://osf.io/3f2h6/
Code
df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |> 
  select(-COUNTRY_OTHER) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |> 
  select(-COUNTRY_OTHER, -EDUCATION_OTHER, -Sample) |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  
  
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |> 
  select(-ID) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
    Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
    Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
    Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
    Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
  )
  • Lin (2023)
    • https://osf.io/3eztd
    • Note: No tickle because same Chinese character
Code
df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) 


df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH)
  • VonMohr (2023) - Study 3
    • https://osf.io/7p9u5/
Code
df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
  mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
  rename(
    Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  • Makowski
    • Note: First sample has some missing items
    • Note: Analog scales
    • note df10c = primals data
Code
df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
  ) |>
  filter(!is.na(Urinate))

df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    typeSample = Sample
  )

df10c <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionPrimals/refs/heads/main/data/data_participants.csv") |>
  select("Participant" = "participant_id", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |>
  filter(!is.na(Heart)) # participant is outlier and had total IAS scores below 0.4 
Code
pi_vars <- c("PI_Enticing", "PI_Alive", "PI_Safe", "PI_Good", "PI_Changing", 
             "PI_Hierarchical", "PI_Understandable", "PI_AboutMe", "PI_Abundant", 
             "PI_Acceptable", "PI_Beautiful", "PI_Cooperative", "PI_Funny", 
             "PI_Harmless", "PI_Improvable", "PI_Intentional", "PI_Interconnected", 
             "PI_Interesting", "PI_Just", "PI_Meaningful", "PI_NeedsMe", 
             "PI_Pleasurable", "PI_Progressing", "PI_Regenerative", "PI_Stable", 
             "PI_WorthExploring")


correlation::correlation(
  df10c[, "IAS_Total", drop = FALSE],  
  data2 = select(df10c, all_of(pi_vars)), 
) |> 
  # Formatting correlation results
  mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
  
  # Releveling factors for visualization
  mutate(Parameter2 = fct_relevel(Parameter2, rev(pi_vars)),
         Parameter1 = fct_relevel(Parameter1, "IAS_Total")) |>
  
  # Plotting the correlation matrix
  ggplot(aes(x = Parameter1, y = Parameter2)) +
  geom_tile(aes(fill = r), color = "white") +  
  geom_text(aes(label = val), size = 3) +  
  labs(title = "Correlation Matrix") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", limit = c(-1, 1), midpoint = 0, guide = guide_colourbar(ticks = FALSE)) +  
  theme_minimal() +
  theme(
    legend.title = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(hjust = 1)
  )

  • Poreiro et al.,
    • Study 1: https://osf.io/49wbv
    • Study 2: unpublished
Code
df11a <- haven::read_sav("../data/Giulia/Interoceptive Attention ESM.sav") |> 
  select(AGE, GENDER, contains("IAS_ACC"), -IAS_ACC_20, starts_with(c("MAIA", "TAS20", "DEP", "ANX")),  -ANX_21, -TAS20_21) |> #anx_21 , tas20_21 are attention checks
  mutate_all(as.numeric) |>
  rename(
    Age = AGE, Gender = GENDER,
    Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4,
    Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7, Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9,
    Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
    Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18,
    Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_21, Itch = IAS_ACC_22
  ) |>
   mutate(Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    Gender %in% c(3, 5) ~ "Other",
    Gender == 4 ~ "NA"
  )) |> ## compute dimensions for the TAS
  mutate(TAS_DIF = TAS20_1 + TAS20_3 + TAS20_6 + TAS20_7 + TAS20_9 + TAS20_13 + TAS20_14,
         TAS_DDF = TAS20_2 + (1-TAS20_4) + TAS20_11 + TAS20_12 + TAS20_17,
         TAS_EOT = (1-TAS20_5) + TAS20_8 + (1-TAS20_10) + TAS20_15 + TAS20_16 + (1-TAS20_18) + (1-TAS20_19) + TAS20_20)


df11b <- haven::read_sav("../data/Giulia/Sleep and Exteroceptive Interoceptive Sensitivity.sav") |> # 165 participants
  filter(Total_Attention_Fail == 0) |> # 32 participants removed based on failing one or more checks
  select(AGE, GENDER, contains("IAS_ACC"), ) |> 
  mutate_all(as.numeric) |> 
  rename( Age = AGE, Gender = GENDER, 
          Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4, Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7,
          Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9, Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
          Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18, Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_20, 
          Itch = IAS_ACC_21 )  |> 
  mutate(Gender = case_when(
    Gender == 1 ~ "Male", 
    Gender == 2 ~"Female", 
    Gender %in% c(3, 5) ~ "Other", 
    Gender == 4 ~ "NA", )) |>
  na.omit()


# Save all
data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b, df10c=df10c, df11a=df11a, df11b=df11b)
save(data, file = "../data/data.RData")

Participants

  • Sample 1a: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
  • Sample 1b: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
  • Sample 2: Data from Gaggero’s(2020) study, downloaded from OSF, included 814 participants (Mean age = 24.9, SD = 5.3, range: [18, 58], 0.2% missing; Gender: 60.3% women, 39.4% men, 0.25% non-binary).
  • Sample 3: Data from Campos’s(2022) study, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
  • Sample 4: Data from Todd’s(2022) study, downloaded from OSF, included 802 participants ().
  • Sample 5: Data from Arslanova (2022) study, downloaded from OSF, included 143 participants (Mean age = 28.5, SD = 7.6, range: [18, 73]; Gender: 45.5% women, 54.5% men, 0.00% non-binary).
  • Sample 6: Data from Brand’s(2022) study, downloaded from OSF, included 619 participants (Mean age = 43.9, SD = 14.5, range: [18, 78]; Gender: 78.7% women, 20.2% men, 1.13% non-binary).
  • Sample 7a: Data from Brand’s(2023) study, downloaded from OSF, included 522 participants (Mean age = 23.4, SD = 6.7, range: [18, 79]; Gender: 79.5% women, 19.7% men, 0.77% non-binary).
  • Sample 7b: Data from Brand’s(2023) study, downloaded from OSF, included 1993 participants (Mean age = 32.0, SD = 12.6, range: [16, 81]; Gender: 77.7% women, 21.7% men, 0.60% non-binary).
  • Sample 7c: Data from Brand’s(2023) study, downloaded from OSF, included 808 participants (Mean age = 27.3, SD = 9.3, range: [18, 72], 0.5% missing; Gender: 68.9% women, 30.2% men, 0.87% non-binary).
  • Sample 8a: Data from Lin’s(2023) study, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
  • Sample 8b: Data from Lin’s(2023) study, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
  • Sample 9: Data from VonMohr’s(2023) study 3, downloaded from OSF, included 21843 participants (Mean age = 56.5, SD = 14.4, range: [18, 93], 0.2% missing; Gender: 73.2% women, 25.1% men, 1.55% non-binary, 0.15% missing).
  • Sample 10a: Data from [Makowski’s(2023)] study , downloaded from OSF, included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).
  • Sample 10b: Data from [Makowski’s(2023)] study , downloaded from OSF, included 836 participants (Mean age = 25.1, SD = 11.3, range: [18, 76], 0.1% missing; Gender: 73.8% women, 22.6% men, 2.87% non-binary, 0.72% missing; Education: Bachelor, 22.85%; Doctorate, 2.15%; High School, 63.52%; Master, 6.22%; missing, 0.24%; Other, 5.02%).
  • Sample 10c: Data from [Makowski’s(2023)] study , downloaded from OSF, included 104 participants (Mean age = 21.6, SD = 5.0, range: [18, 50], 1.0% missing; Gender: 76.0% women, 19.2% men, 2.88% non-binary, 1.92% missing).
  • Sample 11a: Data from Giulia’s (2024) study , included 107 participants (Mean age = 26.8, SD = 9.2, range: [18, 57]; Gender: 74.8% women, 23.4% men, 1.87% non-binary)
  • Sample 11b: Data from [Giulia] study , included 131 participants (Mean age = 30.9, SD = 12.0, range: [18, 60]; Gender: 76.3% women, 22.9% men, 0.76% non-binary)

Total N = 32214.

Code
desc <- function(df){
  summary <- df |> summarise(
      n = n(),
      mean = mean(Age, na.rm = TRUE),
      sd = sd(Age, na.rm = TRUE),
      min = min(Age, na.rm = TRUE),
      max = max(Age, na.rm = TRUE),
      female = sum(Gender == "Female", na.rm = TRUE)
    )
    return(summary)
}

desc_df1a <- desc(df1a)
desc_df1b <- desc(df1b)
desc_df2 <- desc(df2)
desc_df3 <- desc(df3)
#desc_df4 <- desc(df4)
desc_df5 <- desc(df5_attention)
desc_df6 <- desc(df6)
desc_df7a <- desc(df7a)
desc_df7b <- desc(df7b)
desc_df7c <- desc(df7c)
desc_df8a <- desc(df8a)
desc_df8b <- desc(df8b)
desc_df9 <- desc(df9)
desc_df10a <- desc(df10a)
desc_df10b <- desc(df10b)
desc_df10c <- desc(df10c)
desc_df11a <- desc(df11a)
desc_df11b <- desc(df11b)

desc_total <- rbind(desc_df1a, desc_df1b, desc_df2, desc_df3, desc_df5, desc_df6, desc_df7a, desc_df7b, desc_df7c, desc_df8a, desc_df8b, desc_df9, desc_df10a, desc_df10b, desc_df10c, desc_df11a, desc_df11b)

## calculate weighted mean

total_mean = sum(desc_total$mean * desc_total$n) / sum(desc_total$n)
total_sd = sum(desc_total$sd * desc_total$n) / sum(desc_total$n)
perce_female = sum(desc_total$female)/sum(desc_total$n) *100
Code
library(gt)

# APA style ####
gt_apastyle <- function(gt_table, font.size=12) {
  gt_table  |> 
    gt::opt_table_lines(extent = "none") %>%
    gt::tab_options(
      heading.border.bottom.width = 2,
      heading.border.bottom.color = "black",
      heading.border.bottom.style = "solid",
      table.border.top.color = "black",
      table.border.top.style = "solid",
      table.border.top.width = 2,  
      table_body.hlines.color = "white",
      table_body.border.top.color = "black",
      table_body.border.top.style = "solid",
      table_body.border.top.width = 2,
      heading.title.font.size = font.size,
      table.font.size = font.size,
      heading.subtitle.font.size = font.size,
      table_body.border.bottom.color = "black",
      table_body.border.bottom.width = 2,
      table_body.border.bottom.style = "solid",
      column_labels.border.bottom.color = "black",
      column_labels.border.bottom.style = "solid",
      column_labels.border.bottom.width = 1
    ) |> 
      gt::opt_table_font(font = "times")
}

make_age <- function(age) {
  age <- as.numeric(age) 
  
  mean(age, na.rm = TRUE) |> 
    insight::format_value(digits=1) |> 
    paste0(" ± ", 
           insight::format_value(sd(age, na.rm = TRUE), digits=1)) 
} 

# Create the dataframe
table <- data.frame(
   Sample = c('Murphy et al., (2020)', '', '', 'Gaggero et al., (2021)', 'Campos et al., (2022)', 'Todd et al., (2022)', 'Arslanova et al., (2022)', 'Brand et al., (2022)', 'Brand et al., (2023)', '', '', '', 'Lin et al., (2023)', '', '', 'VonMohr et al., (2023)', 'Makowski et al., (2023a)', 'Makowski et al., (2023b)', 'Makowski et al., (2023c)', 'Poreiro et al., (2024)', 'Poreiro et al., unpublished', 'Total'),
  Subsample = c('', 'Sample 1', 'Sample 2', '', '', '', '', '', '', 'Sample 1', 'Sample 2', 'Sample 3', '', 'Sample 1', 'Sample 2', '', '', '', '', '','', ''),   
  Language = c('', 'English', 'English' , 'English and Italian', 'Portuguese', 'English', 'English', 'German', '', 'German', 'German', 'German', '', 'Chinese', 'Chinese', 'English', 'English', 'English', 'English', 'English', 'English',''),
  N = c('', nrow(df1a), nrow(df1b), nrow(df2), nrow(df3), nrow(df4), nrow(df5_attention), nrow(df6),'', nrow(df7a), nrow(df7b), 802,'', nrow(df8a), nrow(df8b), nrow(df9), nrow(df10a), nrow(df10b), nrow(df10c), nrow(df11a), nrow(df11b), 32216),
  Age = c('', make_age(df1a$Age), make_age(df1b$Age), make_age(df2$Age), make_age(df3$Age), "48.6.6 ± 14.1*", make_age(df5_attention$Age), make_age(df6$Age),'', make_age(df7a$Age), make_age(df7b$Age), make_age(df7c$Age),'', make_age(df8a$Age), make_age(df8b$Age), make_age(df9$Age), make_age(df10a$Age), make_age(df10b$Age), make_age(df10c$Age), make_age(df11a$Age), make_age(df11b$Age), '48.6 ± 13.1'),
  Range = c('', '18-69', '18-91', '18-58', '18-72', '18-92*', '18-73', '18-78', '', '18-79', '16-81', '18-72','', '16-60', '20-60', '18-93', '18-73', '17-76','18-50', '18-57', '18-60', '17-93'),
  Female_Percentage = c('', '69.4%', '70.1%', '60.3%', '59.6%', '50%*', '46.8%', '78.7%', '', '79.5%', '77.7%', '68.9%', '', '57.0%', '56.2%', '73.2%', '50.3%', '53.0%', '76%', '74.8%', '75.9%', '71.6'),
  Availability= c('', '', 'https://osf.io/3m5nh', 'https://osf.io/5x9sg', 'https://osf.io/j6ef3', 'https://osf.io/ms354', 'https://osf.io/mp3cy', 'https://osf.io/xwz6g', '', '', 'https://osf.io/3f2h6', '','','', 'https://osf.io/3eztd', 'https://osf.io/7p9u5', 'https://github.com/RealityBending/IllusionGameReliability', 'https://github.com/DominiqueMakowski/PHQ4R', 'https://github.com/RealityBending/InteroceptionPrimals', 'https://osf.io/49wbv', '','')
) 

table_apa <- table |> 
  gt() |>
  cols_align(align = c("right"), columns = "Age") |> 
  cols_label(Age = "Age (Mean  ± SD)", Female_Percentage = "Female %") |> 
  tab_footnote("* Information taken from the sample description of relevant paper rather than recomputed.") |> 
  gt_apastyle() 
  
table_apa
gtsave(table_apa, "figures/table1.tex")

Descriptive

Distribution

The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.

Code
vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 10a",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10b",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10c <- estimate_density(select(df10c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10c",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens11a <- estimate_density(select(df11a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11a")
dens11b <- estimate_density(select(df11b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11b")


p1 <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b, dens10c, dens11a, dens11b) |>
  mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b", "Sample 10c", "Sample 11a", "Sample 11b")) |>
  # mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line() +
  scale_color_material_d() +
  scale_linetype_manual(values = c("solid", "dashed"), guide="none") +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        plot.title = element_text(face="bold")) +
  labs(title = "Distribution of responses for all items across various datasets", x = "Response", y = "Density", color = "Item") +
  facet_wrap(~Sample, scales = "free_y", nrow = 5)

ggsave("figures/Figure1.png", p1, width=10, height=10, dpi=200, bg="white")
p1

Code
data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))
data10c <- select(df10c, all_of(dens10c$Parameter))
data11a <- normalize(select(df11a, all_of(dens11a$Parameter)), range = c(1, 5))
data11b <- na.omit(normalize(select(df11b, all_of(dens11b$Parameter)), range = c(1, 5))) # remove NAS
row.names(data11b) <- NULL



data_all <- rbind(
  data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
  mutate(data8a, Tickle = NA), mutate(data8b, Tickle = NA), data9,
  mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA),  data10b,  data10c, data11a, data11b
) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

Code
make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Parameter1, y = Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
}

make_correlation(data_all)

Code
make_correlation(data1a)

Code
make_correlation(data1b)

Code
make_correlation(data2)

Code
make_correlation(data3)

Code
make_correlation(data4)

Code
make_correlation(data5)

Code
make_correlation(data6)

Code
make_correlation(data7a)

Code
make_correlation(data7b)

Code
make_correlation(data7c)

Code
make_correlation(data8a)

Code
make_correlation(data8b)

Code
make_correlation(data9)

Code
make_correlation(data10a)

Code
make_correlation(data10b)

Code
make_correlation(data10c)

Code
make_correlation(data11a)

Code
make_correlation(data11b)

Graph Analysis

Exploratory Graph Analysis (EGA) is a recently developed framework for psychometric assessment, that can be used to estimate the number of dimensions in questionnaire data using network estimation methods and community detection algorithms, and assess the stability of dimensions and items using bootstrapping. See https://r-ega.net/articles/ega.html for details.

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

Code
uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.364

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.290
 Urinate Defecate 0.265

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Sneeze     Cough 0.230
  Heart Breathing 0.225
 Hungry   Thirsty 0.217
Code
uva$keep_remove
$keep
[1] "Itch"

$remove
[1] "Tickle"
Code
uva <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"
Code
uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208
Code
uva$keep_remove
$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"
Code
uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.415
 Tickle      Itch 0.351

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.293
 Urinate Defecate 0.282
    Wind     Burp 0.275
  Sneeze    Cough 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i node_j   wto
 Muscles   Pain 0.205
Code
uva$keep_remove
$keep
[1] "Heart" "Itch" 

$remove
[1] "Breathing" "Tickle"   
Code
uva <- EGAnet::UVA(data = data5, cut.off = 0.3)
Warning: Some variables did not belong to a dimension: Blood_Sugar, Affective_touch 

Use caution: These variables have been removed from the TEFI calculation
Code
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

    node_i node_j   wto
    Tickle   Itch 0.247
    Bruise   Itch 0.246
 Breathing  Vomit 0.207
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237
Code
uva$keep_remove
$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  
Code
uva <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.434
 Tickle   Itch 0.376

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.295
 Urinate Defecate 0.294

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
$keep
[1] "Sneeze" "Tickle"

$remove
[1] "Cough" "Itch" 
Code
uva <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.317
  Tickle     Itch 0.311

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.281
 Sneeze     Cough 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.237
Code
uva$keep_remove
$keep
[1] "Defecate" "Itch"    

$remove
[1] "Urinate" "Tickle" 
Code
uva <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.457
 Sneeze  Cough 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.280
 Urinate Defecate 0.254

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Hungry   Thirsty 0.242
  Vomit    Sneeze 0.224
  Heart Breathing 0.218
Code
uva$keep_remove
$keep
[1] "Cough" "Itch" 

$remove
[1] "Sneeze" "Tickle"
Code
uva <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.379
    Wind     Burp 0.373
 Urinate Defecate 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.285
 Sneeze     Cough 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j  wto
 Hungry Thirsty 0.24
Code
uva$keep_remove
$keep
[1] "Defecate" "Wind"     "Itch"    

$remove
[1] "Urinate" "Burp"    "Tickle" 
Code
uva <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data10c, cut.off = 0.3)
Warning: Some variables did not belong to a dimension: Heart, Vomit, Blood_Sugar 

Use caution: These variables have been removed from the TEFI calculation
Code
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data11a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.330
   Heart Breathing 0.314
    Wind      Burp 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i  node_j   wto
 Muscles    Pain 0.279
  Tickle    Itch 0.277
  Hungry Thirsty 0.255

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
  Vomit Sneeze 0.201
Code
uva$keep_remove
$keep
[1] "Heart"    "Defecate" "Burp"    

$remove
[1] "Breathing" "Urinate"   "Wind"     
Code
uva <- EGAnet::UVA(data = data11b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Hungry Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j   wto
   Wind    Burp 0.214
   Temp Muscles 0.202
Code
uva$keep_remove
NULL

Networks

Code
make_egas <- function(data) {
  egas <- list()
  for (model in c("glasso")) {  # , "TMFG"
    for (algo in c("walktrap", "louvain")) {
      for (type in c("ega", "ega.fit")) { # "" # "hierega"
        if (type == "ega.fit" & algo == "louvain") next # Too slow
        egas[[paste0(model, "_", algo, "_", type)]] <- EGAnet::bootEGA(
          data = data,
          seed = 123,
          model = model,
          algorithm = algo,
          EGA.type = type,
          type = "resampling",
          plot.itemStability = FALSE,
          verbose = FALSE
        )
      }
    }
  }

  p <- EGAnet::compare.EGA.plots(
    egas$glasso_walktrap_ega, 
    # egas$glasso_walktrap_ega.fit,
    egas$glasso_louvain_ega, 
    # egas$TMFG_louvain_ega,
    # egas$TMFG_walktrap_ega, 
    # egas$TMFG_walktrap_ega.fit,
    labels = c(
      "glasso_walktrap_ega", 
      # "glasso_walktrap_ega.fit",
      "glasso_louvain_ega" 
      # "TMFG_louvain_ega",
      # "TMFG_walktrap_ega", 
      # "TMFG_walktrap_ega.fit"
    ),
    rows = 3,
    plot.all = FALSE
  )$all
  list(egas = egas, p = p)
}


egas_all <- make_egas(data_all)
egas_all$p

Code
egas1a <- make_egas(data1a)
egas1a$p

Code
egas1b <- make_egas(data1b)
egas1b$p

Code
egas2 <- make_egas(data2)
egas2$p

Code
egas3 <- make_egas(data3)
egas3$p

Code
egas4 <- make_egas(data4)
egas4$p

Code
egas5 <- make_egas(data5)
egas5$p

Code
egas6 <- make_egas(data6)
egas6$p

Code
egas7a <- make_egas(data7a)
egas7a$p

Code
egas7b <- make_egas(data7b)
egas7b$p

Code
egas7c <- make_egas(data7c)
egas7c$p

Code
egas8a <- make_egas(data8a)
egas8a$p

Code
egas8b <- make_egas(data8b)
egas8b$p

Code
egas9 <- make_egas(data9)
egas9$p

Code
egas10a <- make_egas(data10a)
egas10a$p

Code
egas10b <- make_egas(data10b)
egas10b$p

Code
egas10c <- make_egas(data10c)
Warning: Some variables did not belong to a dimension: Heart, Vomit, Blood_Sugar 

Use caution: These variables have been removed from the TEFI calculation
Warning in order(as.numeric(names(dimension_stability))): NAs introduced by
coercion
Warning in order(as.numeric(names(dimension_stability))): NAs introduced by
coercion
Code
egas10c$p

Code
egas11a <- make_egas(data11a)
egas11a$p

Code
egas11b <- make_egas(data11b)
egas11b$p

Structure Stability

Code
patchwork::wrap_plots(lapply(egas_all$egas, plot), nrow = 3)
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

Code
patchwork::wrap_plots(lapply(egas1a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas1b$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas2$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas3$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas4$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas5$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas6$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas7a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas7b$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas7c$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas8a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas8b$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas9$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas10a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas10b$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas10c$egas, plot), nrow = 3)
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Code
patchwork::wrap_plots(lapply(egas11a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas11b$egas, plot), nrow = 3)

Scores

Code
ega <- egas_all$egas$glasso_louvain_ega$EGA
plot(ega)

Code
get_scores <- function(data, ega) {
  EGAnet::net.scores(data, ega)$scores$std.scores |>
    as.data.frame() |>
    setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))
}

ega_scores_1a <- get_scores(data1a, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_1b <- get_scores(data1b, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_2 <- get_scores(data2, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_3 <- get_scores(data3, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_4 <- get_scores(data4, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_5 <- get_scores(data5, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_6 <- get_scores(data6, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_7a <- get_scores(data7a, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_7b <- get_scores(data7b, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_7c <- get_scores(data7c, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_8a <- get_scores(data8a, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_8b <- get_scores(data8b, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_9 <- get_scores(data9, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_10a <- get_scores(data10a, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_10b <- get_scores(data10b, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_10c <- get_scores(data10c, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_11a <- get_scores(data11a, ega)
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
ega_scores_11b <- na.omit(get_scores(data11b, ega)) # remove two rows with NAS
The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`
Code
row.names(ega_scores_11b) <- NULL

Factor Analysis

How Many Factors

Code
n <- parameters::n_factors(data_all, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data1a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
3 BIC BIC
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data1b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
3 CNG CNG
3 Scree (SE) Scree_SE
3 BIC BIC
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
5 BIC (adjusted) BIC
Code
n <- parameters::n_factors(data2, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
4 beta Multiple_regression
4 Scree (SE) Scree_SE
4 BIC BIC
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
Code
n <- parameters::n_factors(data3, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
3 CNG CNG
3 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
4 BIC BIC
Code
n <- parameters::n_factors(data4, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
5 BIC BIC
Code
n <- parameters::n_factors(data5, n_max = 12)
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Code
plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
1 BIC BIC
3 CNG CNG
3 BIC (adjusted) BIC
4 beta Multiple_regression
4 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data6, n_max = 12)
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Code
plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data7a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
3 BIC BIC
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data7b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
4 beta Multiple_regression
4 Scree (SE) Scree_SE
5 Bentler Bentler
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
Code
n <- parameters::n_factors(data7c, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data8a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
4 beta Multiple_regression
4 Scree (SE) Scree_SE
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
5 BIC BIC
Code
n <- parameters::n_factors(data8b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
4 beta Multiple_regression
4 BIC BIC
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
5 Scree (SE) Scree_SE
5 BIC (adjusted) BIC
Code
n <- parameters::n_factors(data9, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data10a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 BIC BIC
4 BIC (adjusted) BIC
Code
n <- parameters::n_factors(data10b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
3 Scree (SE) Scree_SE
3 BIC BIC
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data10c, n_max = 12)
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Code
plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
1 BIC BIC
3 Bartlett Barlett
3 CNG CNG
4 beta Multiple_regression
5 Optimal coordinates Scree
5 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data11a, n_max = 12)
Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
The estimated weights for the factor scores are probably incorrect.  Try a
different factor score estimation method.
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
The estimated weights for the factor scores are probably incorrect.  Try a
different factor score estimation method.
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Code
plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
1 BIC BIC
3 Bentler Bentler
3 CNG CNG
3 VSS complexity 2 VSS
4 beta Multiple_regression
5 Optimal coordinates Scree
Code
n <- parameters::n_factors(data11b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
3 CNG CNG
3 Scree (SE) Scree_SE
3 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree

EFA

Code
efa5 <- parameters::factor_analysis(data_all, n = 4, rotation = "oblimin", sort = TRUE)
Loading required namespace: GPArotation
Code
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Burp 0.74 0.03 -0.02 -0.07 1.02 0.48
Cough 0.69 -0.02 0.03 0.06 1.02 0.48
Wind 0.67 0.02 -0.02 1.92e-03 1.00 0.54
Sneeze 0.63 -0.05 0.05 0.14 1.13 0.52
Vomit 0.43 0.06 0.04 0.19 1.43 0.63
Temp 0.26 0.25 0.05 0.21 3.01 0.61
Sex_arousal 0.25 0.23 0.02 0.20 2.94 0.67
Taste 0.24 0.20 0.10 0.14 2.99 0.71
Breathing 0.05 0.61 -0.04 0.03 1.03 0.59
Hungry -0.11 0.54 -5.65e-03 0.20 1.35 0.63
Heart 0.08 0.49 0.03 -0.07 1.11 0.73
Thirsty -0.07 0.47 -5.43e-03 0.25 1.59 0.64
Pain 0.15 0.40 0.10 0.10 1.55 0.62
Muscles 0.30 0.37 0.07 0.02 2.01 0.59
Blood_Sugar 0.11 0.35 0.21 -0.18 2.42 0.76
Bruise 0.24 0.34 0.22 -0.19 3.26 0.65
Tickle -0.04 -0.04 0.82 0.04 1.02 0.38
Itch 0.02 5.03e-03 0.75 -0.04 1.01 0.43
Affective_touch 0.07 0.20 0.33 0.11 2.04 0.70
Urinate 0.04 0.10 0.04 0.65 1.06 0.45
Defecate 0.16 0.04 0.05 0.60 1.15 0.47

The 4 latent factors (oblimin rotation) accounted for 41.58% of the total variance of the original data (MR1 = 14.45%, MR3 = 11.76%, MR2 = 8.09%, MR4 = 7.28%).

Code
efa4 <- parameters::factor_analysis(data1a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR2 MR3 MR1 MR4 Complexity Uniqueness
Tickle 0.72 0.03 -0.10 0.11 1.09 0.47
Itch 0.67 0.14 -9.44e-03 -0.10 1.13 0.46
Affective_touch 0.49 -0.13 0.29 0.13 1.97 0.59
Blood_Sugar 0.29 0.21 0.18 -0.23 3.51 0.75
Burp -0.03 0.73 -0.03 -1.21e-03 1.00 0.50
Cough 7.56e-03 0.69 7.00e-03 -0.01 1.00 0.52
Sneeze 0.07 0.54 7.09e-03 0.20 1.30 0.55
Wind 0.08 0.46 0.07 0.16 1.37 0.62
Hungry -0.16 0.04 0.65 0.03 1.13 0.60
Thirsty -4.84e-03 0.01 0.49 0.11 1.11 0.70
Breathing 0.14 0.05 0.44 0.03 1.24 0.70
Heart 0.11 -0.07 0.40 0.08 1.32 0.78
Bruise 0.31 0.25 0.34 -0.25 3.74 0.59
Pain 0.29 0.01 0.32 0.20 2.67 0.61
Temp 0.19 0.13 0.29 0.16 2.91 0.67
Muscles 0.20 0.13 0.26 0.16 3.24 0.69
Sex_arousal 0.18 0.02 0.21 0.15 2.81 0.82
Defecate 0.10 0.11 2.77e-03 0.66 1.10 0.46
Urinate -0.11 0.13 0.24 0.52 1.65 0.54
Taste 0.23 0.06 0.17 0.30 2.60 0.69
Vomit 0.21 0.18 0.11 0.25 3.26 0.71

The 4 latent factors (oblimin rotation) accounted for 37.94% of the total variance of the original data (MR2 = 10.46%, MR3 = 10.17%, MR1 = 10.08%, MR4 = 7.23%).

Code
efa4 <- parameters::factor_analysis(data1b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR3 MR4 Complexity Uniqueness
Urinate 0.69 -0.18 0.13 -0.04 1.22 0.54
Defecate 0.68 -0.13 0.04 0.06 1.10 0.55
Breathing 0.56 0.03 0.06 7.96e-03 1.03 0.63
Hungry 0.55 0.15 -7.15e-03 0.02 1.16 0.60
Thirsty 0.52 0.05 0.14 -0.01 1.16 0.62
Sex_arousal 0.50 0.14 -0.04 0.20 1.49 0.56
Temp 0.46 0.23 0.16 -0.01 1.77 0.52
Pain 0.46 0.43 -0.05 -0.03 2.03 0.49
Taste 0.40 0.21 -0.07 0.12 1.77 0.69
Heart 0.40 0.06 -0.04 0.15 1.36 0.76
Muscles 0.35 0.24 0.07 0.09 2.07 0.65
Vomit 0.32 0.06 0.23 0.11 2.20 0.68
Tickle 0.02 0.69 0.01 0.04 1.01 0.47
Itch -0.02 0.67 0.18 0.03 1.16 0.41
Bruise 0.07 0.44 0.15 0.06 1.32 0.65
Blood_Sugar -0.08 0.38 0.16 0.09 1.61 0.77
Affective_touch 0.38 0.38 -0.13 0.05 2.26 0.63
Sneeze 0.05 0.05 0.71 0.01 1.02 0.41
Cough 0.04 0.02 0.71 0.04 1.01 0.43
Wind 7.52e-03 -0.04 -5.18e-03 0.94 1.00 0.15
Burp -0.01 0.27 0.17 0.47 1.91 0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Code
efa4 <- parameters::factor_analysis(data2, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR4 MR2 Complexity Uniqueness
Urinate 0.62 0.12 0.03 -0.16 1.21 0.55
Sex_arousal 0.58 0.02 0.03 -3.44e-04 1.01 0.64
Hungry 0.50 -0.07 0.02 0.25 1.52 0.65
Defecate 0.50 0.25 0.04 -0.12 1.62 0.56
Thirsty 0.47 -6.03e-03 0.07 0.14 1.22 0.70
Breathing 0.38 6.09e-03 -0.02 0.31 1.93 0.71
Temp 0.35 0.26 0.06 0.11 2.13 0.63
Vomit 0.33 0.22 0.06 0.05 1.84 0.72
Heart 0.32 -0.08 0.02 0.28 2.10 0.80
Taste 0.28 0.28 0.10 -0.06 2.32 0.73
Cough -7.64e-03 0.69 0.04 0.09 1.04 0.46
Sneeze 0.15 0.58 0.02 -0.06 1.16 0.57
Burp 0.12 0.42 0.09 0.15 1.51 0.63
Wind 0.19 0.42 0.07 0.08 1.52 0.63
Tickle 0.02 -0.01 0.83 -0.06 1.01 0.35
Itch -0.11 0.23 0.51 0.19 1.80 0.54
Affective_touch 0.24 -0.19 0.37 0.19 2.89 0.72
Bruise -0.06 0.23 0.05 0.53 1.40 0.58
Blood_Sugar -0.04 -7.08e-03 0.22 0.40 1.56 0.75
Pain 0.20 0.23 0.02 0.36 2.33 0.63
Muscles 0.19 0.21 -0.01 0.27 2.75 0.75

The 4 latent factors (oblimin rotation) accounted for 36.65% of the total variance of the original data (MR1 = 12.95%, MR3 = 10.42%, MR4 = 6.75%, MR2 = 6.52%).

Code
efa4 <- parameters::factor_analysis(data3, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Cough 0.88 -0.08 -0.04 0.11 1.05 0.29
Sneeze 0.74 0.05 1.42e-03 -0.03 1.01 0.41
Vomit 0.55 0.12 0.04 -0.09 1.17 0.58
Burp 0.55 0.08 0.22 -0.13 1.47 0.48
Wind 0.53 0.10 0.21 -0.25 1.84 0.48
Temp 0.46 0.34 -0.06 0.09 1.97 0.52
Taste 0.29 0.18 0.10 0.12 2.37 0.74
Hungry -8.85e-03 0.70 -0.01 -0.09 1.03 0.54
Urinate 0.08 0.66 -0.12 0.11 1.15 0.52
Thirsty -0.08 0.64 -0.03 0.14 1.14 0.62
Defecate 0.10 0.63 -3.50e-05 -0.07 1.08 0.53
Breathing -0.10 0.49 0.28 -0.07 1.76 0.66
Sex_arousal 0.17 0.44 0.17 -0.22 2.20 0.58
Pain 0.21 0.41 0.08 0.18 2.00 0.58
Heart 0.05 0.28 0.24 0.08 2.23 0.77
Muscles 0.24 0.27 0.20 0.19 3.66 0.61
Tickle 0.02 -0.05 0.84 -0.04 1.01 0.31
Itch 0.03 0.02 0.74 0.11 1.05 0.39
Affective_touch 0.08 0.23 0.36 0.12 2.08 0.67
Bruise 0.17 0.07 0.29 0.45 2.08 0.55
Blood_Sugar 4.11e-03 0.11 0.35 0.43 2.07 0.60

The 4 latent factors (oblimin rotation) accounted for 45.69% of the total variance of the original data (MR1 = 15.50%, MR3 = 15.46%, MR2 = 11.02%, MR4 = 3.71%).

Code
efa4 <- parameters::factor_analysis(data4, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR4 MR2 MR3 Complexity Uniqueness
Burp 0.69 -0.09 -4.33e-03 0.04 1.04 0.57
Cough 0.66 0.03 0.04 -1.08e-03 1.01 0.52
Wind 0.64 1.45e-03 0.03 0.02 1.01 0.56
Sneeze 0.52 0.07 0.07 -3.65e-03 1.07 0.64
Vomit 0.50 0.05 0.02 0.10 1.10 0.66
Pain 0.40 0.28 0.04 0.12 2.00 0.56
Muscles 0.39 0.23 0.12 0.07 1.91 0.60
Sex_arousal 0.29 0.28 0.08 0.02 2.18 0.70
Bruise 0.24 -0.13 0.20 0.23 3.47 0.77
Hungry -0.12 0.62 0.04 0.17 1.25 0.56
Urinate 0.20 0.60 -0.02 -0.09 1.27 0.55
Thirsty -0.10 0.53 0.09 0.15 1.30 0.66
Defecate 0.26 0.48 -0.04 -0.01 1.58 0.61
Temp 0.27 0.40 0.07 0.04 1.84 0.60
Taste 0.20 0.36 0.15 0.04 1.99 0.65
Itch 5.98e-03 -0.07 0.82 0.01 1.01 0.35
Tickle 8.91e-04 0.06 0.73 -0.05 1.03 0.47
Affective_touch 0.04 0.24 0.45 0.04 1.58 0.63
Blood_Sugar 0.06 -0.13 0.33 0.19 2.05 0.81
Breathing 2.50e-03 0.09 -0.03 0.78 1.03 0.34
Heart 0.06 -0.08 0.03 0.70 1.05 0.50

The 4 latent factors (oblimin rotation) accounted for 41.38% of the total variance of the original data (MR1 = 14.47%, MR4 = 10.83%, MR2 = 8.95%, MR3 = 7.14%).

Code
efa4 <- parameters::factor_analysis(data5, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR4 MR2 MR3 Complexity Uniqueness
Sex_arousal 0.53 0.02 -0.05 -0.10 1.09 0.73
Wind 0.52 0.14 0.06 0.06 1.19 0.60
Heart 0.50 -0.08 0.04 0.08 1.12 0.75
Taste 0.44 -0.03 -0.02 0.17 1.29 0.76
Muscles 0.42 -0.03 0.10 0.15 1.38 0.75
Hungry 0.38 0.06 -0.17 -3.75e-03 1.42 0.84
Burp 0.38 0.06 0.14 -0.17 1.78 0.80
Affective_touch 0.36 -4.99e-03 -0.07 0.03 1.09 0.87
Defecate 0.35 0.23 -0.14 0.03 2.08 0.76
Thirsty 0.34 0.22 -0.03 -0.10 1.98 0.78
Blood_Sugar 0.27 -0.18 0.05 -0.01 1.85 0.93
Sneeze 0.27 0.17 0.26 -0.08 2.91 0.75
Urinate -0.05 0.72 -0.11 0.06 1.07 0.51
Cough 0.04 0.63 0.19 0.01 1.19 0.49
Pain 0.12 0.57 0.10 -0.10 1.21 0.59
Temp 0.29 0.30 -0.06 0.09 2.27 0.72
Itch -0.04 0.01 0.94 0.03 1.01 0.12
Bruise 0.20 0.13 0.43 -0.03 1.64 0.68
Tickle 0.27 -0.18 0.38 0.04 2.30 0.75
Breathing 3.71e-03 3.31e-03 0.02 0.96 1.00 0.07
Vomit 0.25 0.26 -4.07e-03 0.29 2.97 0.67

The 4 latent factors (oblimin rotation) accounted for 33.72% of the total variance of the original data (MR1 = 12.07%, MR4 = 8.86%, MR2 = 7.16%, MR3 = 5.62%).

Code
efa4 <- parameters::factor_analysis(data6, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR4 MR2 MR3 Complexity Uniqueness
Burp 0.79 0.03 -0.05 0.02 1.01 0.38
Cough 0.66 0.06 0.17 -0.28 1.53 0.32
Wind 0.62 0.10 0.03 0.03 1.06 0.49
Bruise 0.58 0.04 -0.13 0.16 1.25 0.66
Breathing 0.55 -0.05 0.09 0.22 1.38 0.59
Muscles 0.49 0.20 0.06 0.09 1.42 0.53
Sneeze 0.49 0.08 0.32 -0.26 2.37 0.37
Blood_Sugar 0.49 0.02 -0.15 0.18 1.49 0.76
Heart 0.43 -0.04 7.61e-03 0.24 1.59 0.75
Vomit 0.43 0.09 0.16 -0.15 1.64 0.63
Temp 0.34 0.20 0.27 -0.09 2.72 0.52
Taste 0.25 0.20 0.18 0.10 3.13 0.69
Tickle -0.05 0.89 -0.06 0.02 1.02 0.32
Itch 0.02 0.74 0.08 -0.10 1.06 0.37
Affective_touch 0.14 0.48 0.01 0.14 1.36 0.62
Pain 0.24 0.43 0.07 0.12 1.81 0.54
Sex_arousal 0.24 0.29 0.14 0.14 2.93 0.62
Defecate 0.04 0.02 0.79 5.46e-03 1.01 0.32
Urinate -0.06 0.05 0.75 0.06 1.03 0.44
Hungry 0.13 0.12 0.28 0.47 1.94 0.54
Thirsty 0.10 0.03 0.38 0.44 2.09 0.58

The 4 latent factors (oblimin rotation) accounted for 47.47% of the total variance of the original data (MR1 = 20.28%, MR4 = 12.37%, MR2 = 10.55%, MR3 = 4.27%).

Code
efa4 <- parameters::factor_analysis(data7a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Itch 0.63 -0.12 0.22 0.07 1.34 0.46
Pain 0.59 -0.01 0.03 0.09 1.05 0.59
Tickle 0.52 -0.12 0.14 0.14 1.41 0.61
Muscles 0.49 0.13 -0.03 0.07 1.19 0.67
Affective_touch 0.45 0.13 0.08 -2.36e-03 1.24 0.70
Breathing 0.41 0.26 -0.16 -0.03 2.03 0.76
Bruise 0.38 0.35 -0.04 -0.05 2.05 0.68
Temp 0.33 -0.04 0.17 0.20 2.22 0.71
Blood_Sugar 0.33 0.33 0.02 -0.13 2.30 0.75
Heart 0.30 0.30 -0.06 -0.22 2.88 0.83
Burp -0.04 0.66 0.20 0.07 1.22 0.39
Wind -0.01 0.61 0.08 0.15 1.17 0.50
Cough 0.03 0.11 0.76 -7.60e-04 1.05 0.31
Sneeze 0.07 5.91e-03 0.70 8.40e-05 1.02 0.47
Vomit 0.11 0.25 0.25 0.18 3.18 0.67
Urinate 3.91e-04 0.06 2.65e-04 0.57 1.02 0.65
Defecate 0.06 0.07 0.04 0.56 1.06 0.61
Hungry 0.30 5.52e-03 -0.23 0.36 2.69 0.72
Sex_arousal 0.13 0.23 1.59e-03 0.34 2.09 0.71
Thirsty 0.24 0.12 -0.10 0.27 2.69 0.79
Taste 0.17 0.12 0.19 0.21 3.52 0.75

The 4 latent factors (oblimin rotation) accounted for 36.52% of the total variance of the original data (MR1 = 13.12%, MR3 = 8.41%, MR2 = 7.80%, MR4 = 7.18%).

Code
efa4 <- parameters::factor_analysis(data7b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Sneeze 0.68 0.07 0.04 -0.16 1.14 0.53
Cough 0.67 0.06 0.04 -0.07 1.04 0.51
Burp 0.65 -0.07 -0.02 0.19 1.21 0.50
Wind 0.53 -0.01 2.01e-03 0.20 1.28 0.59
Vomit 0.48 0.05 0.05 -7.54e-03 1.05 0.72
Temp 0.25 0.23 0.15 0.09 2.91 0.67
Muscles 0.25 0.24 0.06 0.20 3.05 0.65
Sex_arousal 0.24 0.08 0.21 0.16 3.04 0.70
Tickle -0.06 0.75 -0.02 -3.98e-03 1.01 0.50
Itch 0.05 0.72 0.02 -0.02 1.01 0.43
Affective_touch 0.07 0.42 0.01 0.14 1.31 0.70
Pain 0.13 0.34 0.09 0.21 2.15 0.63
Urinate -0.04 5.49e-03 0.75 -6.11e-03 1.01 0.46
Defecate 0.06 -1.14e-03 0.69 -0.04 1.02 0.50
Thirsty -5.07e-03 -0.04 0.35 0.24 1.83 0.79
Taste 0.16 0.13 0.26 0.17 3.07 0.69
Breathing 0.02 0.11 0.09 0.47 1.19 0.68
Heart 0.03 0.12 -0.01 0.45 1.16 0.74
Blood_Sugar 0.08 0.05 0.02 0.38 1.14 0.79
Bruise 0.23 0.10 8.46e-03 0.34 1.98 0.70
Hungry -3.08e-03 0.02 0.28 0.32 1.97 0.75

The 4 latent factors (oblimin rotation) accounted for 36.96% of the total variance of the original data (MR1 = 12.41%, MR3 = 9.28%, MR2 = 8.10%, MR4 = 7.16%).

Code
efa4 <- parameters::factor_analysis(data7c, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR3 MR1 MR4 MR2 Complexity Uniqueness
Cough 0.73 -0.04 0.04 0.04 1.02 0.44
Sneeze 0.73 -4.69e-03 -0.03 9.49e-03 1.00 0.48
Wind 0.56 0.10 0.09 -2.95e-03 1.12 0.58
Vomit 0.50 0.04 -0.01 -6.16e-03 1.02 0.74
Burp 0.47 0.13 0.12 0.08 1.34 0.60
Urinate 0.10 0.56 -0.05 4.96e-03 1.09 0.65
Defecate 0.24 0.55 -0.11 -0.03 1.47 0.60
Hungry -0.16 0.54 0.14 0.03 1.31 0.69
Sex_arousal 0.09 0.45 0.08 0.04 1.15 0.71
Thirsty -0.06 0.39 0.20 0.03 1.56 0.76
Temp 0.18 0.27 0.21 0.05 2.82 0.71
Taste 0.12 0.25 0.22 0.05 2.53 0.77
Bruise 0.11 -0.09 0.54 0.12 1.26 0.61
Muscles 0.13 0.06 0.46 -0.10 1.32 0.72
Breathing 0.04 0.22 0.42 -0.08 1.62 0.71
Blood_Sugar 0.02 -0.12 0.41 0.10 1.30 0.81
Pain 0.02 0.29 0.39 0.06 1.90 0.64
Heart 0.08 0.07 0.31 0.05 1.29 0.84
Tickle -0.01 0.06 -0.06 0.85 1.02 0.30
Itch 0.05 -0.08 0.07 0.78 1.04 0.34
Affective_touch -0.05 0.29 0.12 0.31 2.37 0.75

The 4 latent factors (oblimin rotation) accounted for 35.87% of the total variance of the original data (MR3 = 11.03%, MR1 = 9.41%, MR4 = 7.85%, MR2 = 7.57%).

Code
efa4 <- parameters::factor_analysis(data8a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Sneeze 0.69 -0.03 -0.01 0.03 1.01 0.53
Cough 0.68 -0.07 0.03 0.04 1.04 0.54
Burp 0.55 0.18 -2.31e-04 -0.09 1.26 0.61
Wind 0.48 0.08 0.07 -0.03 1.09 0.70
Vomit 0.42 0.09 -0.04 0.13 1.31 0.74
Taste 0.20 0.19 0.18 0.04 3.08 0.79
Bruise -0.03 0.65 -9.02e-03 0.04 1.01 0.58
Affective_touch 0.07 0.55 -0.03 -0.02 1.05 0.67
Blood_Sugar 0.15 0.45 -0.09 0.09 1.40 0.70
Itch 0.07 0.41 0.14 0.04 1.30 0.71
Muscles 0.02 0.39 0.15 0.02 1.29 0.77
Thirsty 5.49e-03 -0.05 0.67 -0.02 1.01 0.58
Hungry 0.06 -0.10 0.49 0.17 1.35 0.67
Pain -0.04 0.31 0.41 0.02 1.89 0.67
Urinate 0.13 0.12 0.39 -2.22e-03 1.43 0.73
Temp 0.03 0.16 0.36 0.09 1.58 0.75
Defecate 0.24 0.09 0.34 -0.04 2.00 0.72
Heart -0.02 0.04 -0.04 0.65 1.01 0.58
Breathing 0.08 -6.25e-03 0.12 0.49 1.17 0.66
Sex_arousal 0.06 0.19 0.17 0.22 3.04 0.78

The 4 latent factors (oblimin rotation) accounted for 32.61% of the total variance of the original data (MR1 = 10.53%, MR3 = 9.06%, MR2 = 8.24%, MR4 = 4.78%).

Code
efa4 <- parameters::factor_analysis(data8b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR4 MR2 MR3 Complexity Uniqueness
Burp 0.65 0.07 6.06e-03 -0.10 1.07 0.56
Cough 0.63 -0.06 0.03 0.11 1.09 0.56
Sneeze 0.57 0.02 -8.72e-03 0.13 1.10 0.61
Wind 0.54 0.02 -0.03 1.35e-04 1.01 0.71
Defecate 0.34 0.11 0.28 -0.13 2.48 0.70
Vomit 0.27 0.18 0.06 0.13 2.39 0.77
Bruise -0.06 0.64 -0.04 0.01 1.02 0.63
Affective_touch 0.08 0.55 -0.07 -0.08 1.11 0.67
Blood_Sugar 0.09 0.50 -0.03 0.02 1.07 0.71
Muscles -2.84e-04 0.41 0.12 0.05 1.19 0.79
Itch 0.15 0.40 0.02 0.13 1.52 0.69
Pain 0.05 0.29 0.26 0.02 2.07 0.78
Taste 0.09 0.24 0.21 0.11 2.76 0.79
Thirsty 0.03 -0.10 0.71 -0.04 1.05 0.51
Hungry -4.08e-03 0.04 0.51 0.17 1.23 0.65
Temp -0.09 0.22 0.39 0.10 1.87 0.77
Urinate 0.19 0.18 0.33 -0.11 2.49 0.73
Breathing 0.06 -0.03 -0.01 0.70 1.02 0.50
Heart -0.01 0.07 0.04 0.50 1.06 0.72
Sex_arousal 0.06 0.17 0.17 0.30 2.35 0.76

The 4 latent factors (oblimin rotation) accounted for 32.02% of the total variance of the original data (MR1 = 9.96%, MR4 = 9.20%, MR2 = 7.36%, MR3 = 5.49%).

Code
efa4 <- parameters::factor_analysis(data9, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Burp 0.78 -0.01 -9.54e-03 -0.06 1.01 0.45
Wind 0.71 -0.03 -7.51e-03 0.03 1.01 0.51
Cough 0.69 0.02 0.01 0.04 1.01 0.47
Sneeze 0.63 -1.87e-03 0.02 0.11 1.07 0.51
Vomit 0.42 0.06 0.06 0.19 1.49 0.63
Temp 0.30 0.24 0.03 0.20 2.75 0.60
Sex_arousal 0.27 0.21 0.03 0.18 2.76 0.68
Taste 0.24 0.19 0.12 0.13 3.07 0.72
Breathing 0.04 0.62 -0.03 0.02 1.02 0.58
Hungry -0.09 0.59 -0.01 0.15 1.17 0.60
Thirsty -0.06 0.51 -0.01 0.23 1.42 0.60
Heart 0.08 0.47 0.06 -0.08 1.16 0.74
Pain 0.18 0.41 0.09 0.06 1.52 0.62
Muscles 0.34 0.37 0.06 7.80e-04 2.05 0.57
Blood_Sugar 0.08 0.31 0.28 -0.16 2.68 0.76
Bruise 0.24 0.29 0.28 -0.19 3.70 0.66
Tickle -0.04 -0.04 0.79 0.05 1.02 0.41
Itch 0.01 -0.01 0.77 -0.02 1.00 0.41
Affective_touch 0.05 0.19 0.34 0.11 1.85 0.71
Urinate 0.02 0.09 0.03 0.69 1.04 0.41
Defecate 0.11 0.02 0.05 0.67 1.06 0.43

The 4 latent factors (oblimin rotation) accounted for 42.44% of the total variance of the original data (MR1 = 14.92%, MR3 = 11.60%, MR2 = 8.42%, MR4 = 7.50%).

Code
efa4 <- parameters::factor_analysis(data10a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR4 MR2 Complexity Uniqueness
Hungry 0.62 -0.04 0.15 0.04 1.13 0.51
Sex_arousal 0.58 0.06 0.01 0.03 1.03 0.60
Pain 0.55 0.04 0.09 -8.09e-03 1.06 0.61
Thirsty 0.55 0.09 0.07 0.03 1.10 0.58
Temp 0.36 0.11 0.19 0.09 1.85 0.65
Affective_touch 0.28 -0.05 0.13 0.26 2.51 0.75
Defecate 0.06 0.67 -7.30e-03 0.02 1.02 0.49
Sneeze -0.18 0.53 0.29 0.06 1.85 0.59
Vomit 0.13 0.53 -0.05 0.09 1.20 0.61
Urinate 0.40 0.47 -0.12 -0.02 2.09 0.52
Burp 2.70e-03 0.46 0.17 0.08 1.34 0.64
Breathing 0.02 0.02 0.71 0.01 1.00 0.45
Heart 0.17 7.44e-03 0.59 0.01 1.16 0.51
Muscles 0.15 0.07 0.47 0.07 1.31 0.60
Itch -5.06e-03 -0.02 -0.02 0.78 1.00 0.43
Tickle -0.01 0.02 -0.01 0.73 1.00 0.47
Bruise 8.63e-03 0.04 0.06 0.38 1.08 0.81
Wind 0.10 0.08 0.10 0.26 1.82 0.82

The 4 latent factors (oblimin rotation) accounted for 40.83% of the total variance of the original data (MR1 = 12.34%, MR3 = 10.09%, MR4 = 9.29%, MR2 = 9.12%).

Code
efa4 <- parameters::factor_analysis(data10b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR4 MR3 MR1 MR2 Complexity Uniqueness
Breathing 0.55 -0.14 0.10 0.14 1.32 0.62
Muscles 0.55 0.16 0.02 -0.02 1.18 0.59
Heart 0.48 -0.01 0.04 0.08 1.07 0.71
Taste 0.42 0.20 -0.05 0.09 1.55 0.67
Temp 0.38 0.05 0.38 -0.06 2.09 0.56
Pain 0.36 0.07 0.32 -0.04 2.07 0.63
Cough 0.30 0.27 -0.07 0.16 2.62 0.69
Defecate 0.07 0.60 0.07 -0.06 1.07 0.58
Wind -0.03 0.57 0.05 0.17 1.19 0.57
Burp 0.12 0.39 -0.08 0.27 2.13 0.64
Vomit 0.03 0.31 0.18 0.11 1.90 0.77
Sneeze 0.19 0.31 -0.06 0.23 2.64 0.70
Hungry 4.71e-03 0.04 0.63 0.07 1.03 0.55
Thirsty 0.21 -0.02 0.50 0.11 1.44 0.55
Urinate -0.03 0.35 0.42 -0.03 1.96 0.62
Affective_touch 0.08 0.07 0.39 0.18 1.59 0.69
Sex_arousal 0.17 0.32 0.34 -0.10 2.64 0.61
Itch 0.15 0.09 -0.10 0.54 1.30 0.59
Bruise -0.07 -0.03 0.23 0.53 1.39 0.67
Tickle 0.06 0.15 0.06 0.48 1.26 0.63
Blood_Sugar 0.05 -0.12 0.20 0.32 2.07 0.84

The 4 latent factors (oblimin rotation) accounted for 35.81% of the total variance of the original data (MR4 = 10.18%, MR3 = 9.19%, MR1 = 9.17%, MR2 = 7.27%).

Code
efa4 <- parameters::factor_analysis(data10c, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR4 MR3 Complexity Uniqueness
Hungry 0.67 -0.11 -0.15 0.20 1.34 0.53
Thirsty 0.63 0.08 -0.04 -0.05 1.06 0.59
Sex_arousal 0.57 0.14 0.14 -0.07 1.28 0.53
Pain 0.54 0.07 0.14 0.07 1.20 0.57
Urinate 0.49 0.04 0.35 -0.08 1.89 0.52
Temp 0.35 0.04 0.35 0.15 2.37 0.56
Muscles 0.33 0.10 0.33 -0.09 2.35 0.66
Affective_touch 0.28 0.18 -0.02 0.13 2.15 0.81
Itch -0.05 0.79 -0.01 -0.06 1.02 0.41
Bruise 0.07 0.55 -0.07 0.18 1.29 0.63
Taste 0.28 0.42 -0.14 0.09 2.10 0.67
Tickle 0.07 0.40 0.28 -0.06 1.90 0.66
Breathing 0.20 0.25 0.22 0.09 3.16 0.72
Wind -0.05 -0.04 0.69 0.09 1.05 0.53
Sneeze 0.29 0.02 0.35 0.23 2.72 0.57
Defecate 0.25 0.20 0.28 0.04 2.85 0.68
Heart 0.04 0.12 0.22 0.04 1.66 0.90
Vomit 0.15 -0.01 0.17 0.04 2.11 0.92
Burp 0.03 -0.03 0.03 0.75 1.01 0.41
Blood_Sugar -0.25 0.20 0.17 0.35 2.97 0.82
Cough 0.05 0.26 0.07 0.34 2.00 0.75

The 4 latent factors (oblimin rotation) accounted for 35.93% of the total variance of the original data (MR1 = 13.46%, MR2 = 8.82%, MR4 = 8.11%, MR3 = 5.54%).

Code
efa4 <- parameters::factor_analysis(data11a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR4 MR3 Complexity Uniqueness
Cough 0.78 -0.02 -1.87e-03 0.15 1.07 0.34
Heart 0.57 -0.10 -0.04 -0.22 1.36 0.69
Breathing 0.56 0.14 0.07 -0.22 1.50 0.59
Sneeze 0.44 0.16 0.07 0.19 1.73 0.64
Wind 0.40 0.40 0.02 0.14 2.25 0.54
Temp 0.35 -0.08 0.32 -0.08 2.21 0.76
Vomit 0.32 0.11 0.17 -0.05 1.83 0.80
Defecate 8.51e-03 0.76 0.08 -0.09 1.05 0.38
Urinate -0.08 0.67 0.19 -5.85e-03 1.19 0.45
Burp 0.35 0.45 -0.16 0.27 2.91 0.53
Sex_arousal 0.37 0.37 0.10 -0.10 2.28 0.60
Muscles 0.27 0.28 0.11 0.16 2.97 0.70
Thirsty -0.13 0.16 0.74 0.08 1.19 0.36
Taste 0.14 0.12 0.49 -0.03 1.31 0.64
Blood_Sugar 0.30 -0.19 0.47 6.77e-03 2.09 0.68
Hungry 0.06 0.13 0.45 -0.07 1.25 0.72
Affective_touch 0.31 -0.23 0.36 0.08 2.82 0.74
Pain 0.11 0.29 0.33 0.12 2.51 0.65
Itch -0.04 -0.04 -5.51e-03 0.72 1.01 0.49
Tickle 0.05 -0.03 0.12 0.52 1.13 0.69
Bruise 0.06 -0.17 0.18 0.41 1.79 0.78

The 4 latent factors (oblimin rotation) accounted for 39.22% of the total variance of the original data (MR1 = 12.83%, MR2 = 10.45%, MR4 = 9.59%, MR3 = 6.35%).

Code
efa4 <- parameters::factor_analysis(data11b, n = 4, rotation = "oblimin", sort = TRUE)
Warning in GPFoblq(A, Tmat = Tmat, normalize = normalize, eps = eps, maxit =
maxit, : convergence not obtained in GPFoblq. 1000 iterations used.
Code
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR2 MR1 MR4 MR3 Complexity Uniqueness
Hungry 0.75 -0.09 0.03 0.06 1.05 0.46
Thirsty 0.74 0.07 -0.19 0.06 1.16 0.43
Breathing 0.71 0.08 -0.19 0.14 1.26 0.42
Urinate 0.67 0.01 0.16 -0.21 1.31 0.47
Pain 0.65 6.54e-03 0.35 0.10 1.61 0.28
Defecate 0.51 0.18 0.20 -0.28 2.22 0.47
Muscles 0.41 0.20 0.04 0.31 2.40 0.52
Temp 0.38 0.38 0.19 -0.02 2.45 0.43
Heart 0.37 0.12 0.19 0.06 1.79 0.70
Burp 0.08 0.83 -0.20 0.20 1.27 0.26
Sneeze 0.05 0.69 0.21 -0.32 1.64 0.26
Cough 0.02 0.66 0.09 -0.01 1.04 0.49
Wind 0.11 0.61 -0.02 0.06 1.10 0.55
Itch -0.22 0.59 0.35 0.10 2.01 0.45
Vomit 0.12 0.51 0.22 -0.13 1.61 0.50
Tickle 0.10 0.26 0.50 0.05 1.62 0.50
Taste 0.14 0.19 0.47 -0.07 1.60 0.57
Affective_touch 0.10 0.12 0.41 0.35 2.29 0.60
Sex_arousal 0.22 0.17 0.41 0.14 2.24 0.57
Blood_Sugar 0.05 0.02 0.16 0.55 1.18 0.67
Bruise 0.27 0.17 0.11 0.50 1.92 0.48

The 4 latent factors (oblimin rotation) accounted for 51.95% of the total variance of the original data (MR2 = 19.08%, MR1 = 17.78%, MR4 = 9.34%, MR3 = 5.75%).

CFA

We discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant.

Ambiguous: - Temp - Vomit - Affective_touch - Sex_arousal - Taste

Lower-Level Factors

Code
m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"

m4 <- "
Sustenance =~ Hungry + Thirsty + Urinate + Defecate + Muscles + Pain
Skin =~ Itch + Bruise
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m5 <- "
Sustenance =~ Hungry + Thirsty + Pain
UrinateDefecate =~ Urinate + Defecate
Skin =~ Itch + Bruise + Muscles
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m6 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"

m7 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
CoughSneeze =~ Cough + Sneeze
WindBurp =~ Wind + Burp
"

is_converged <- function(m) {
  ev <- eigen(lavaan::lavTech(m, "cov.lv")[[1]], symmetric=TRUE, only.values=TRUE)$values
  if (any(ev < -1 * .Machine$double.eps^(3/4))) {
    "Not converged"
  } else {
    ""
  }
}

compare_mods <- function(data) {
  
  fit_model <- function(f, data) {
    if(!"Cough" %in% names(data)) {
      f <- str_remove(f, fixed("Cough +"))
    }
    lavaan::cfa(f, data = data)
  }
  
  
  cfa1 <- fit_model(m1, data = data)
  cfa4 <- fit_model(m4, data = data)
  cfa5 <- fit_model(m5, data = data)
  cfa6 <- fit_model(m6, data = data)
  cfa7 <- fit_model(m7, data = data)
  
  anova(cfa1, cfa4, cfa5, cfa6, cfa7) |> 
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(cfa1, cfa4, cfa5, cfa6, cfa7), is_converged))
}

display(compare_mods(data_all))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -171638 -171228 2195.20 56
cfa6 -167039 -166680 6805.73 62 < .001
cfa5 -162764 -162446 11091.21 67 < .001
cfa4 -157684 -157400 16178.96 71 < .001
cfa1 -144334 -144100 29540.73 77 < .001
Code
m7_all <- lavaan::cfa(m7, data = data_all)
parameters::parameters(m7_all, standardize=TRUE) |>
  display()
# Loading
Link Coefficient SE 95% CI z p
HungryThirsty =~ Hungry 0.68 4.69e-03 (0.67, 0.69) 145.39 < .001
HungryThirsty =~ Thirsty 0.71 4.68e-03 (0.70, 0.72) 151.38 < .001
UrinateDefecate =~ Urinate 0.77 3.88e-03 (0.76, 0.78) 198.08 < .001
UrinateDefecate =~ Defecate 0.77 3.88e-03 (0.76, 0.78) 199.00 < .001
ItchBruise =~ Itch 0.50 5.65e-03 (0.49, 0.52) 89.31 < .001
ItchBruise =~ Bruise 0.67 5.89e-03 (0.66, 0.68) 113.21 < .001
MusclesPain =~ Muscles 0.71 4.10e-03 (0.70, 0.72) 173.39 < .001
MusclesPain =~ Pain 0.68 4.17e-03 (0.67, 0.69) 163.56 < .001
HeartBreathing =~ Breathing 0.79 5.12e-03 (0.78, 0.80) 153.73 < .001
HeartBreathing =~ Heart 0.59 5.06e-03 (0.58, 0.60) 116.02 < .001
CoughSneeze =~ Cough 0.80 3.47e-03 (0.80, 0.81) 230.97 < .001
CoughSneeze =~ Sneeze 0.75 3.62e-03 (0.74, 0.75) 206.40 < .001
WindBurp =~ Wind 0.77 3.55e-03 (0.76, 0.77) 216.49 < .001
WindBurp =~ Burp 0.82 3.42e-03 (0.81, 0.82) 239.00 < .001
# Correlation
Link Coefficient SE 95% CI z p
HungryThirsty ~~ UrinateDefecate 0.65 6.15e-03 (0.64, 0.66) 105.82 < .001
HungryThirsty ~~ ItchBruise 0.45 8.89e-03 (0.43, 0.47) 50.88 < .001
HungryThirsty ~~ MusclesPain 0.63 6.91e-03 (0.62, 0.65) 91.58 < .001
HungryThirsty ~~ HeartBreathing 0.64 6.90e-03 (0.63, 0.66) 93.27 < .001
HungryThirsty ~~ CoughSneeze 0.46 7.03e-03 (0.45, 0.47) 65.42 < .001
HungryThirsty ~~ WindBurp 0.43 7.02e-03 (0.42, 0.45) 61.80 < .001
UrinateDefecate ~~ ItchBruise 0.44 8.27e-03 (0.42, 0.46) 53.09 < .001
UrinateDefecate ~~ MusclesPain 0.65 6.14e-03 (0.64, 0.66) 105.65 < .001
UrinateDefecate ~~ HeartBreathing 0.53 6.71e-03 (0.52, 0.55) 79.61 < .001
UrinateDefecate ~~ CoughSneeze 0.56 5.95e-03 (0.55, 0.57) 94.13 < .001
UrinateDefecate ~~ WindBurp 0.51 6.10e-03 (0.50, 0.52) 83.40 < .001
ItchBruise ~~ MusclesPain 0.80 7.86e-03 (0.78, 0.81) 101.69 < .001
ItchBruise ~~ HeartBreathing 0.51 8.62e-03 (0.49, 0.53) 59.31 < .001
ItchBruise ~~ CoughSneeze 0.64 7.60e-03 (0.63, 0.65) 84.23 < .001
ItchBruise ~~ WindBurp 0.64 7.46e-03 (0.63, 0.66) 86.06 < .001
MusclesPain ~~ HeartBreathing 0.64 6.88e-03 (0.63, 0.65) 93.18 < .001
MusclesPain ~~ CoughSneeze 0.67 6.00e-03 (0.65, 0.68) 110.95 < .001
MusclesPain ~~ WindBurp 0.63 6.04e-03 (0.62, 0.65) 105.01 < .001
HeartBreathing ~~ CoughSneeze 0.52 6.73e-03 (0.50, 0.53) 76.91 < .001
HeartBreathing ~~ WindBurp 0.46 6.86e-03 (0.44, 0.47) 66.38 < .001
CoughSneeze ~~ WindBurp 0.72 4.80e-03 (0.71, 0.73) 150.93 < .001
Code
display(compare_mods(data1a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -1951 -1749 88.22 56
cfa6 -1925 -1748 126.29 62 < .001
cfa5 -1860 -1703 201.38 67 < .001 Not converged
cfa4 -1834 -1694 234.82 71 < .001
cfa1 -1661 -1546 420.24 77 < .001
Code
display(compare_mods(data1b))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -1898 -1706 66.50 56
cfa6 -1823 -1654 154.19 62 < .001
cfa5 -1801 -1652 185.98 67 < .001
cfa4 -1755 -1621 240.22 71 < .001
cfa1 -1632 -1523 374.42 77 < .001
Code
display(compare_mods(data2))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -4588 -4358 123.77 56
cfa6 -4565 -4363 158.82 62 < .001
cfa5 -4490 -4311 243.72 67 < .001
cfa4 -4426 -4266 315.85 71 < .001
cfa1 -4259 -4127 495.03 77 < .001
Code
display(compare_mods(data3))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -3588 -3380 114.03 56
cfa6 -3508 -3325 206.45 62 < .001
cfa5 -3453 -3291 271.41 67 < .001
cfa4 -3419 -3275 313.05 71 < .001
cfa1 -3143 -3024 600.82 77 < .001
Code
display(compare_mods(data4))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -3397 -3167 81.35 56
cfa6 -3321 -3120 168.76 62 < .001
cfa5 -3142 -2964 357.91 67 < .001
cfa4 -3069 -2910 438.71 71 < .001
cfa1 -2627 -2496 892.62 77 < .001
Code
display(compare_mods(data5))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -96 49 101.12 56
cfa6 -101 26 107.73 62 0.359
cfa5 -102 10 116.58 67 0.115 Not converged
cfa4 -104 -3 123.39 71 0.146 Not converged
cfa1 -84 -1 154.66 77 < .001 Not converged
Code
display(compare_mods(data6))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -4359 -4142 152.88 56
cfa6 -4146 -3956 376.94 62 < .001 Not converged
cfa5 -4061 -3893 472.34 67 < .001 Not converged
cfa4 -3865 -3714 676.52 71 < .001 Not converged
cfa1 -3699 -3575 854.37 77 < .001 Not converged
Code
display(compare_mods(data7a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -3268 -3059 106.90 56
cfa6 -3133 -2950 254.03 62 < .001 Not converged
cfa5 -3115 -2953 282.30 67 < .001 Not converged
cfa4 -3053 -2909 351.61 71 < .001 Not converged
cfa1 -2842 -2723 574.94 77 < .001 Not converged
Code
display(compare_mods(data7b))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -11892 -11617 188.21 56
cfa6 -11703 -11462 388.89 62 < .001 Not converged
cfa5 -11583 -11370 518.70 67 < .001
cfa4 -11256 -11066 853.27 71 < .001 Not converged
cfa1 -10705 -10549 1416.33 77 < .001 Not converged
Code
display(compare_mods(data7c))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -1513 -1295 122.56 56
cfa6 -1429 -1237 219.24 62 < .001
cfa5 -1359 -1190 298.55 67 < .001
cfa4 -1296 -1144 370.09 71 < .001
cfa1 -1060 -935 617.75 77 < .001
Code
display(compare_mods(data8a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -5137 -4889 115.89 56
cfa6 -5103 -4885 161.81 62 < .001
cfa5 -5041 -4849 233.76 67 < .001
cfa4 -4930 -4758 352.62 71 < .001
cfa1 -4591 -4449 703.81 77 < .001
Code
display(compare_mods(data8b))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -1738 -1531 83.61 56
cfa6 -1711 -1530 122.14 62 < .001
cfa5 -1683 -1523 160.21 67 < .001
cfa4 -1630 -1487 221.58 71 < .001
cfa1 -1486 -1368 377.91 77 < .001
Code
display(compare_mods(data9))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -137262 -136870 2003.31 56
cfa6 -133150 -132807 6126.47 62 < .001
cfa5 -129698 -129394 9589.14 67 < .001
cfa4 -125505 -125233 13790.36 71 < .001
cfa1 -114370 -114146 24936.87 77 < .001
Code
display(compare_mods(data10a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -1824 -1632 101.11 45
cfa6 -1826 -1654 109.52 50 0.135
cfa5 -1822 -1671 123.48 55 0.016
cfa4 -1752 -1619 200.78 59 < .001
cfa1 -1665 -1556 300.09 65 < .001
Code
display(compare_mods(data10b))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -3965 -3734 185.78 56
cfa6 -3959 -3755 204.27 62 0.005
cfa5 -3938 -3759 234.70 67 < .001 Not converged
cfa4 -3875 -3714 306.23 71 < .001
cfa1 -3710 -3578 482.94 77 < .001
Code
display(compare_mods(data10c))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -432 -303 67.23 56
cfa6 -438 -324 73.36 62 0.409
cfa5 -443 -343 78.32 67 0.421 Not converged
cfa4 -448 -358 81.66 71 0.502 Not converged
cfa1 -434 -360 107.51 77 < .001 Not converged
Code
display(compare_mods(data11a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -219 -88 62.63 56
cfa6 -208 -93 86.39 62 < .001
cfa5 -179 -78 124.90 67 < .001 Not converged
cfa4 -171 -80 141.41 71 0.002
cfa1 -99 -24 224.81 77 < .001
Code
display(compare_mods(data11b))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
cfa7 -612 -471 133.75 56
cfa6 -597 -474 159.97 62 < .001 Not converged
cfa5 -589 -480 178.10 67 0.003 Not converged
cfa4 -587 -489 188.60 71 0.033 Not converged
cfa1 -506 -425 281.32 77 < .001 Not converged

Conclusion: The 7th factor solution was preferred in most datasets, including by indices that penalize increased number of parameters (such as BIC).

Higher Order Factors

Code
m7h1 <- paste0(m7, "
              Interoception =~ HungryThirsty + UrinateDefecate + ItchBruise + MusclesPain + HeartBreathing + CoughSneeze + WindBurp
              ")

m7h2 <- paste0(m7, "
              Expulsion =~ CoughSneeze + WindBurp
              Interoception =~ HungryThirsty + UrinateDefecate + HeartBreathing
              ")

m7h3 <- paste0(m7, "
              Valenced =~ MusclesPain + HeartBreathing + ItchBruise
              Expulsion =~ CoughSneeze + WindBurp
              Homeostasis =~ HungryThirsty + UrinateDefecate
              ")

compare_mods_h <- function(data) {
  
  fit_model_h <- function(f, data) {
    if(!"Cough" %in% names(data)) {
      f <- str_remove(f, fixed("Cough +"))
    }
    lavaan::cfa(f, data = data)
  }
  
  m7 <- fit_model_h(m7, data = data)
  m7h2 <- fit_model_h(m7h2, data = data)
  m7h1 <- fit_model_h(m7h1, data = data)
  m7h3 <- fit_model_h(m7h3, data = data)
  
  anova(m7, m7h1, m7h2, m7h3) |>
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(m7, m7h1, m7h2, m7h3), is_converged))
}

display(compare_mods_h(data_all))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -171638 -171228 2195.20 56
m7h2 -170853 -170527 2999.77 66 < .001
m7h3 -169568 -169250 4286.71 67 < .001
m7h1 -167097 -166805 6763.53 70 < .001

No evidence for higher order factors.

Code
display(compare_mods_h(data1a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -1951 -1749 88.22 56
m7h2 -1930 -1770 128.74 66 < .001
m7h3 -1904 -1747 157.51 67 < .001
m7h1 -1873 -1729 193.71 70 < .001
Code
display(compare_mods_h(data1b))
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -1898 -1706 66.50 56
m7h2 -1901 -1748 83.62 66 0.072
m7h3 -1855 -1706 131.76 67 < .001 Not converged
m7h1 -1841 -1704 151.59 70 < .001
Code
display(compare_mods_h(data2))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -4588 -4358 123.77 56
m7h2 -4573 -4390 158.33 66 < .001
m7h3 -4563 -4385 170.60 67 < .001
m7h1 -4536 -4372 203.47 70 < .001
Code
display(compare_mods_h(data3))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -3588 -3380 114.03 56
m7h2 -3585 -3420 136.81 66 0.012
m7h3 -3581 -3419 143.30 67 0.011
m7h1 -3497 -3348 233.30 70 < .001
Code
display(compare_mods_h(data4))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -3397 -3167 81.35 56
m7h2 -3344 -3161 154.44 66 < .001
m7h3 -3335 -3157 165.03 67 0.001
m7h1 -3282 -3118 223.71 70 < .001
Code
display(compare_mods_h(data5))
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -96 49 101.12 56 Not converged
m7h2 -95 21 122.05 66 0.022 Not converged
m7h3 -95 18 124.10 67 0.152 Not converged
m7h1 -100 3 124.72 70 0.893 Not converged
Code
display(compare_mods_h(data6))
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -4359 -4142 152.88 56 Not converged
m7h2 -4329 -4156 202.72 66 < .001 Not converged
m7h3 -4329 -4161 204.07 67 0.246 Not converged
m7h1 -4311 -4157 227.94 70 < .001 Not converged
Code
display(compare_mods_h(data7a))
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -3268 -3059 106.90 56 Not converged
m7h2 -3260 -3094 134.89 66 0.002 Not converged
m7h3 -3255 -3093 141.68 67 0.009 Not converged
m7h1 -3208 -3059 195.06 70 < .001 Not converged
Code
display(compare_mods_h(data7b))
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -11892 -11617 188.21 56 Not converged
m7h2 -11876 -11658 223.53 66 < .001 Not converged
m7h3 -11872 -11659 229.97 67 0.011 Not converged
m7h1 -11770 -11574 337.40 70 < .001 Not converged
Code
display(compare_mods_h(data7c))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -1513 -1295 122.56 56
m7h2 -1499 -1325 156.56 66 < .001
m7h3 -1483 -1313 175.42 67 < .001
m7h1 -1437 -1281 226.58 70 < .001
Code
display(compare_mods_h(data8a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -5137 -4889 115.89 56
m7h2 -5120 -4923 152.67 66 < .001
m7h3 -5103 -4911 171.44 67 < .001
m7h1 -5007 -4830 273.77 70 < .001
Code
display(compare_mods_h(data8b))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -1738 -1531 83.61 56
m7h2 -1729 -1565 112.49 66 0.001
m7h3 -1727 -1567 116.28 67 0.052
m7h1 -1710 -1562 139.61 70 < .001
Code
display(compare_mods_h(data9))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -137262 -136870 2003.31 56
m7h2 -136603 -136292 2681.71 66 < .001
m7h3 -135349 -135045 3938.35 67 < .001
m7h1 -133209 -132929 6083.88 70 < .001
Code
display(compare_mods_h(data10a))
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -1824 -1632 101.11 45
m7h2 -1798 -1648 146.92 55 < .001
m7h3 -1794 -1647 153.35 56 0.011 Not converged
m7h1 -1778 -1645 174.82 59 < .001 Not converged
Code
display(compare_mods_h(data10b))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -3965 -3734 185.78 56
m7h2 -3938 -3753 233.42 66 < .001
m7h3 -3912 -3732 260.95 67 < .001
m7h1 -3899 -3734 279.76 70 < .001
Code
display(compare_mods_h(data10c))
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -432 -303 67.23 56 Not converged
m7h2 -444 -341 75.22 66 0.630 Not converged
m7h3 -442 -341 79.80 67 0.032 Not converged
m7h1 -444 -351 83.57 70 0.287 Not converged
Code
display(compare_mods_h(data11a))
Warning: lavaan->lav_object_post_check():  
   some estimated ov variances are negative
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   some estimated ov variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated ov variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -219 -88 62.63 56
m7h2 -220 -115 82.41 66 0.031
m7h3 -213 -111 91.40 67 0.003 Not converged
m7h1 -211 -117 99.05 70 0.054
Code
display(compare_mods_h(data11b))
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -612 -471 133.75 56 Not converged
m7h2 -598 -486 167.45 66 < .001
m7h3 -580 -471 187.06 67 < .001 Not converged
m7h1 -555 -454 218.25 70 < .001 Not converged

Model Performance

Code
rbind(
  performance::performance(m7_all, metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "All"),
  performance::performance(update(m7_all, data = data1a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1a"),
  performance::performance(update(m7_all, data = data1b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1b"),
  performance::performance(update(m7_all, data = data2), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 2"),
  performance::performance(update(m7_all, data = data3), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 3"),
  performance::performance(update(m7_all, data = data4), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 4"),
  performance::performance(update(m7_all, data = data5), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 5"),
  performance::performance(update(m7_all, data = data6), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 6"),
  performance::performance(update(m7_all, data = data7a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7a"),
  performance::performance(update(m7_all, data = data7b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7b"),
  performance::performance(update(m7_all, data = data7c), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7c"),
  performance::performance(update(m7_all, data = data8a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8a"),
  performance::performance(update(m7_all, data = data8b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8b"),
  performance::performance(update(m7_all, data = data9), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 9"),
  performance::performance(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data10a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10a"),
  performance::performance(update(m7_all, data = data10b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10b"),
  performance::performance(update(m7_all, data = data10c), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10c"),
  performance::performance(update(m7_all, data = data11a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 11a"),
  performance::performance(update(m7_all, data = data11b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 11b")
) |>
  display()
Chi2 RMSEA CFI SRMR Sample
2195.20 0.03 0.98 0.02 All
88.22 0.04 0.98 0.03 Sample 1a
66.50 0.02 0.99 0.03 Sample 1b
123.77 0.04 0.97 0.03 Sample 2
114.03 0.04 0.98 0.03 Sample 3
81.35 0.02 0.99 0.02 Sample 4
101.12 0.08 0.86 0.07 Sample 5
152.88 0.05 0.97 0.03 Sample 6
106.90 0.04 0.97 0.04 Sample 7a
188.21 0.03 0.98 0.02 Sample 7b
122.56 0.04 0.97 0.03 Sample 7c
115.89 0.03 0.98 0.02 Sample 8a
83.61 0.03 0.98 0.03 Sample 8b
2003.31 0.04 0.98 0.02 Sample 9
82.59 0.05 0.97 0.03 Sample 10a
185.78 0.05 0.95 0.03 Sample 10b
67.23 0.04 0.96 0.06 Sample 10c
62.63 0.03 0.98 0.06 Sample 11a
133.75 0.10 0.90 0.08 Sample 11b

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

Code
# Add empirical variables
add_empirical <- function(data, df) {
  x <- data |>
    mutate(
      HungryThirsty = (df$Hungry + df$Thirsty) / 2,
      UrinateDefecate = (df$Urinate + df$Defecate) / 2,
      MusclesPain = (df$Muscles + df$Pain) / 2,
      HeartBreathing = (df$Heart + df$Breathing) / 2,
      ItchBruise = (df$Itch + df$Bruise) / 2,
      WindBurp = (df$Wind + df$Burp) / 2)
  if("Cough" %in% names(df)) {
    mutate(x, CoughSneeze = (df$Cough + df$Sneeze) / 2)
  } else {
    mutate(x, CoughSneeze = df$Sneeze)
  }
}

# Deal with missing data
temp <- predict(m7_all, newdata=data7c)
scores7c <- data.frame(matrix(ncol = 7, nrow = nrow(ega_scores_7c))) 
scores7c[!is.na(ega_scores_7c$EGA1), ] <- temp
names(scores7c) <- names(as.data.frame(temp))

scores <- list(
  sample1a = cbind(ega_scores_1a, data_addprefix(as.data.frame(predict(m7_all, newdata=data1a)), "CFA_")) |>
    add_empirical(data1a) |> 
    mutate(Sample = "Sample 1a"),
  sample1b = cbind(ega_scores_1b, data_addprefix(as.data.frame(predict(m7_all, newdata=data1b)), "CFA_")) |>
    add_empirical(data1b) |> 
    mutate(Sample = "Sample 1b"),
  sample2 = cbind(ega_scores_2, data_addprefix(as.data.frame(predict(m7_all, newdata=data2)), "CFA_")) |>
    add_empirical(data2) |> 
    mutate(Sample = "Sample 2"),
  sample3 = cbind(ega_scores_3, data_addprefix(as.data.frame(predict(m7_all, newdata=data3)), "CFA_")) |>
    add_empirical(data3) |> 
    mutate(Sample = "Sample 3"),
  sample4 = cbind(ega_scores_4, data_addprefix(as.data.frame(predict(m7_all, newdata=data4)), "CFA_")) |>
    add_empirical(data4) |> 
    mutate(Sample = "Sample 4"),
  sample5 = cbind(ega_scores_5, data_addprefix(as.data.frame(predict(m7_all, newdata=data5)), "CFA_")) |>
    add_empirical(data5) |> 
    mutate(Sample = "Sample 5"),
  sample6 = cbind(ega_scores_6, data_addprefix(as.data.frame(predict(m7_all, newdata=data6)), "CFA_")) |>
    add_empirical(data6) |> 
    mutate(Sample = "Sample 6"),
  sample7a = cbind(ega_scores_7a, data_addprefix(as.data.frame(predict(m7_all, newdata=data7a)), "CFA_")) |>
    add_empirical(data7a) |> 
    mutate(Sample = "Sample 7a"),
  sample7b = cbind(ega_scores_7b, data_addprefix(as.data.frame(predict(m7_all, newdata=data7b)), "CFA_")) |>
    add_empirical(data7b) |> 
    mutate(Sample = "Sample 7b"),
  sample7c = cbind(ega_scores_7c, data_addprefix(scores7c, "CFA_")) |> 
    add_empirical(data7c) |> 
    mutate(Sample = "Sample 7c"),
  sample8a = cbind(ega_scores_8a, data_addprefix(as.data.frame(predict(m7_all, newdata=data8a)), "CFA_")) |>
    add_empirical(data8a) |> 
    mutate(Sample = "Sample 8a"),
  sample8b = cbind(ega_scores_8b, data_addprefix(as.data.frame(predict(m7_all, newdata=data8b)), "CFA_")) |>
    add_empirical(data8b) |> 
    mutate(Sample = "Sample 8b"),
  sample9 = cbind(ega_scores_9, data_addprefix(as.data.frame(predict(m7_all, newdata=data9)), "CFA_")) |>
    add_empirical(data9) |> 
    mutate(Sample = "Sample 9"),
  sample10a = cbind(ega_scores_10a, data_addprefix(as.data.frame(predict(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data_all), newdata=data10a)), "CFA_")) |>
    add_empirical(data10a) |> 
    mutate(Sample = "Sample 10a", CFA_CoughSneeze = NA),
  sample10b = cbind(ega_scores_10b, data_addprefix(as.data.frame(predict(m7_all, newdata=data10b)), "CFA_")) |>
    add_empirical(data10b) |> 
    mutate(Sample = "Sample 10b"),
  sample10c = cbind(ega_scores_10c, data_addprefix(as.data.frame(predict(m7_all, newdata=data10c)), "CFA_")) |>
    add_empirical(data10c) |> 
    mutate(Sample = "Sample 10c"),
  sample11a = cbind(ega_scores_11a, data_addprefix(as.data.frame(predict(m7_all, newdata=data11a)), "CFA_")) |>
    add_empirical(data11a) |> 
    mutate(Sample = "Sample 11a"),
  sample11b = cbind(ega_scores_11b, data_addprefix(as.data.frame(predict(m7_all, newdata=data11b)), "CFA_")) |>
    add_empirical(data11b) |> 
    mutate(Sample = "Sample 11b")
)

save(scores, file = "../data/scores.RData")